Walther etal. BMC Systems Biology 2014, 8:22 
http://www.bionnedcentral.conn/l 752-0509/8/22 



Systems Biology 



SOFTWARE Open Access 



GraTeLPy: graph-theoretic linear stability 
analysis 

Georg R Walther^ ^ Matthew Hartley^ and Maya Mincheva^ 



Abstract 

Background: A biochemical mechanism with mass action l<inetics can be represented as a directed bipartite graph 
(bipartite digraph), and modeled by a system of differential equations. If the differential equations (DE) model can give 
rise to some instability such as multistability or Turing instability, then the bipartite digraph contains a structure 
referred to as a critical fragment. In some cases the existence of a critical fragment indicates that the DE model can 
display oscillations for some parameter values. We have implemented a graph-theoretic method that identifies the 
critical fragments of the bipartite digraph of a biochemical mechanism. 

Results: GraTeLPy lists all critical fragments of the bipartite digraph of a given biochemical mechanism, thus enabling 
a preliminary analysis on the potential of a biochemical mechanism for some instability based on its topological 
structure. The correctness of the implementation is supported by multiple examples. The code is implemented in 
Python, relies on open software, and is available under the GNU General Public License. 

Conclusions: GraTeLPy can be used by researchers to test large biochemical mechanisms with mass action kinetics 
for their capacity for multistability, oscillations and Turing instability. 

Keywords: Biochemical mechanism. Bipartite digraph, Multistability, Turing instability. Oscillations, Parameter-free 
model discrimination 



Background 

Biochemical mechanisms are often modeled by differen- 
tial equations (DE) systems. Instabilities, such as multi- 
stability, oscillations, or Turing instability, are ubiquitous 
in DE models of biochemical mechanisms. Methods from 
bifurcation analysis are usually applied in order to ana- 
lyze DE models for instabilities [1]. Bifurcation analysis 
methods are easily applied when the DE model has one or 
two concentration species (phase plane analysis) or has a 
relatively small number of parameters (numerical bifurca- 
tion analysis). However, it is both difficult and expensive 
to apply bifurcation methods to analyze large DE models 
with many variables for instabilities. 

On the other hand a biochemical mechanism can 
be represented as a directed bipartite graph (bipartite 
digraph), which is a graph with two different sets of nodes 
representing species and reactions, and directed edges 
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connecting a species and a reaction node. The existence 
of structures referred to as critical fragments in the bipar- 
tite digraph of a biochemical mechanism is necessary for 
the existence of multistability or Turing instability in the 
DE model [2-4]. Thus biochemical mechanisms that do 
not have the potential for multistability or Turing insta- 
bility can be ruled out early in the modeling process. 
The existence of a critical fragment that does not contain 
all species nodes can indicate that oscillations exist for 
some parameter values for the DE model [3]. Thus graph- 
theoretic methods can be used to determine the poten- 
tial of various biochemical mechanisms to exhibit some 
desired behavior, including multistability related to cell 
decision [5,6], oscillations related to circadian rhythms 
[7], or Turing instability related to pattern formation [8]. 

Graph-theoretic methods are applicable to mechanisms 
with any number of species and reactions, which enables 
the screening of large biochemical mechanisms for poten- 
tial instabilities. However, application of graph-theoretic 
methods by hand becomes challenging for large mech- 
anisms, making a computational implementation highly 
desirable. 
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The graph-theoretic method implemented by GraTeLPy 
identifies all critical fragments in the bipartite digraph 
of a biochemical mechanism that can give rise to some 
instability (multistability, oscillations and Turing instabil- 
ity) [3,4]. GraTelPy is implemented in Python and can run 
in parallel on computer clusters which increases the size 
of testable biochemical mechanisms. 

Other software packages implement theoretical and 
computational methods for studying chemical reaction 
networks for multistability. Using the deficiency theory 
developed by M. Feinberg and collaborators [9], it can 
be shown that a chemical network model does not admit 
multistability for any choice of parameter values. The 
CRNT toolbox [10] developed originally by M. Feinberg 
implements the Deficiency One algorithm [11], that can 
be used to detect if a given network has the capacity for 
multistability [12]. If a given network admits multiple pos- 
itive equilibria, in many cases the CRNT toolbox returns 
rate constant values such that the corresponding model 
system has at least two positive equilibria. In recent years, 
the CRNT toolbox has been extended to implement an 
algorithm for the mass-action injectivity test. A special 
case of this test is the Jacobian criterion, which provides a 
sufficient condition for excluding the existence of multiple 
positive equilibria and is based on the theory developed in 
[2,13,14]. 

Related software packages include BioNetX [15] and 
CoNtRol [16]. BioNetX is based on the work of M. 
Banaji and G. Craciun [17,18] and is created by C. 
Pantea. BioNetX is used to analyze uni-molecular and 
bi-molecular reaction networks for the existence of 
multiple positive equilibria in [19,20]. CoNtRol [16] 
is a web-based software package that employs matrix 
and graph-theoretic methods based on the DSR graph 
[17,18,21]. In particular CoNtRol provides information 
about the capacity of a given chemical network for multi- 
stability based on the DSR graph and on some additional 
tests. In addition CoNtRol calculates the deficiency of 
a network and checks if a network is weakly reversible. 
BioNetX and CoNtRol are available to download for free, 
they are open-source and are conveniently web-based. 

We describe in Section Mathematical background the 
DE model and the bipartite digraph of a biochemical 
mechanism, as well as the instability criteria. In 
Section Implementation we describe the algorithm 
for finding critical fragments. In Section Results and 
discussion we present several examples along with 
some concluding remarks. A guide for downloading 
and installing GraTeLPy for Mac, Windows and Linux 
operating systems is available in the Additional file 1. 

Mathematical background 

Here we introduce the differential equations model and 
the bipartite digraph representation of a biochemical 



mechanism. In this section we also briefly describe the 
instability criteria for multistability, oscillations and Tur- 
ing instability. More details on the instability criteria are 
available in [3,4]. 

Mathematical model 

A biochemical mechanism withn species A/, / = 1, . . . , fz, 
and m elementary reactions Bj can be written as 

« /c- 

Bj: ; = l,...,m, (1) 

i=l i=l 

where A/ > 0,/ = 1, . . . , m are the rate constants. The con- 
stants aji > 0 and fiji > 0 in (1) are small integers called 
stoichiometric coefficients that account for the number of 
molecules of species At participating in the elementary 
reaction. An example of a biochemical mechanism, the 
reversible substrate inhibition mechanism, is given below: 

Bi: Ai 0, 
B2: 0 ^ Ai, 

^3 : Ai + A2 ^ A3, 
A3 ^ A2, 
^5 : Ai + A3 ^ A4, 
Be: A4^Ai+A3, 

where the first two reactions represent an inflow and 
outflow reaction, respectively. 

We will assume that every species in (1) is consumed 
and produced in at least one true reaction, i.e, a reaction 
which is different from an outflow reaction A]^ 0 or 
an inflow reaction 0 Aj^, However, we do not specif- 
ically require that all species participate in an inflow and 
an outflow reaction. 

We further assume mass action kinetics for the mecha- 
nism (1) with rate functions 

Wj = kjVL^^ . . . M^^", / = 1, . . . , m, (3) 

where is the concentration at time ^ of a species A^, 
/:=!,...,«. 

The ordinary differential equations (ODE) model of a 
mass-action biochemical mechanism (1) can be written in 
vector form as 

u{t) = Sw(u), (4) 

where u{t) = (wi(^), . . . , Un{t))^ is the concentration vec- 
tor of the chemical species of (1), Sji = Pji — aji are 
the entries of the stoichiometric matrix S and w(u) = 
{wi{u)y . . . yWyniu))^ is the vector of rate functions (3). 
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Throughout the paper it will be assumed that the ODE 
system (4) has a positive equilibrium. 

The model equations of the reversible substrate mecha- 
nism (2) are given below 



ui = -kiui + k2 - huiU2 - ksUiU3 + k^u^, 
U2 = -k3UiU2 + k^us, 
Us = k3UiU2 — k^us — ksuius + k^u^, 
U4 = ksuius — keU4, 



(5) 



The rank of the stoichiometric matrix S of (5) equals 3 
since there is one conservation relationship U2 + U3-\-U4 = 
Co. 
Since 

m 

Uiit) =fi(u) = ^SjiWjiu) 

7=1 

where Sji are the stoichiometric matrix entries and wj(u) 
are the rate functions (3), the Jacobian matrix /(u, w) has 
entries that can be written as 



dfi y-^ Wj 

Jikiu.w) = — = ySjiajk — , 
^ ' ' Uk 



(6) 



Note that the concentrations U]^, k = 1, . . . , n and the 
rate functions wj(u), j = 1, . . . , m (both considered eval- 
uated at a positive equilibrium) are used as parameters in 
(6). The rank of the Jacobian (6) equals the rank of the 
stoichiometric matrix S [3]. 

The Jacobian matrix of the model (5) parametrized in 
(u, w) has rank 3 and is given below 



J(u,w) = 



1^1+1^3+1^5 


_W3 


_W5 


we 


Ui 


U2 


U3 




_W3 


_W3 


W4 


0 


Ui 


U2 


Us 


W3-W5 


W3 


14/3 +M/4 




Ui 


U2 


U3 




Ws 


0 


W5 


_we 


Ul 


US 


M4 



(7) 



The characteristic polynomial of J(u, w) is 

n 

P(X) = det(J(u, w) - XI) = Y^akiu, w))J' 

k=o 



(8) 



where / is the identity matrix. Note that the coefficients 
ai = ai(u,w), i = 1, of (8) are also functions of 
(u,w). For example, the last non-zero coefficient of the 
characteristic polynomial of the Jacobian (7) is 



asiu, w) 



UiU^U^, 



U\U2U4 



U\U2U^ 



The bipartite digraph of a biochemical mechanism 

For the convenience of the reader, in this section we 
present definitions regarding the bipartite digraph of a 
biochemical mechanism (1) [3,4,22]. To illustrate the def- 
initions in this section we will continue to use as an 
example the reversible substrate mechanism (2). 

A directed bipartite graph (bipartite digraph) has a node 
set that consists of two disjoint subsets, V\ and V2, and 
each of its directed edges (arcs) has one end in V\ and the 
other in V2 [23]. 

The bipartite digraph G of a biochemical reaction 
network (1) is defined as follows. The nodes are sepa- 
rated into two sets, one for the chemical species V\ = 
{Ai,A2, . . . , A^} and one for the elementary reactions 
V2 = {B\,B2y . . .yB^}, We draw an arc from Aj^ to Bj if 
and only if species Aj^ is a reactant in reaction i.e., if 
the stoichiometric coefficient aj]^ > 0 in (1). Similarly, we 
draw an arc from Bj to Ai if and only if Ai is a product 
in reaction y, i.e., if the stoichiometric coefficient Pji > 0 
in (1). Therefore the set of arcs E{G) consists of arcs such 
as {Aj^yBj) and {Bj.Ai), Hence the bipartite digraph can be 
defined as G = {V,E{G)} where V =ViVJV2\s the set of 
nodes and E{G) is the set of arcs. If an arc is not weighted 
explicitly, we assume that its weight equals 1. The cor- 
responding bipartite digraph of the reversible substrate 
inhibition mechanism (2) is shown in Figure 1. 

The element [Aj^yBj] is an edge if ajk > 0, i.e., if species 
A/^ is a reactant in reaction The weight of an edge E = 
[Ai(,Bj] is defined as 



Ke 



(10) 



(9) 




Figure 1 Bipartite digraph of the reversible substrate inhibition 
mechanism. Bipartite digrapli of tlie reversible reaction meclianism 
(2). Circles denote species nodes and squares denote reaction nodes 
of the mechanism. 
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For example, the edge E =[^1,^3] in Figure 1 has weight 
Ke = -1. 

If aj/^Pji > 0, then the arcs (A/^.Bj) and (Bj^Ai) form 
a positive path [Aj^.Bj.Ai] that corresponds to the pro- 
duction of Ai from Aj^ in 3. reaction The weight of the 
positive path [A^yBj^Ai] is defined as ajj^fiji. For example, 
the positive path [^1,^3, A3] in Figure 1 has weight 1. 

If aji^aji > 0, then the arcs (A/^yBj) and (Ai^Bj) form a 
negative path [A]^,Bj,Ai\ that corresponds to A]^ and Ai 
interacting as reactants in reaction The weight of the 
negative path [A]^,Bj,Ai\ is defined as —ajj^aji. Note that 
the negative paths [Aj^.Bj.Ai] and [Ai,Bj,Aj(] are consid- 
ered to be different since they start at a different species 
node. For example, both [^1,^3,^2] and [^2,^3,^1] in 
Figure 1 are negative paths with weight —1. We note that 
the direction of the arcs is followed in the positive paths 
but not in the negative paths. 

A cycle C of G is a sequence of distinct paths with 
the last species node of each path being the same as the 
first species node of the next path C = {(Ai^^Bj^fAi^), 
{Ai^,Bj^,Ai^),. . ., {Ai^_^,Bjj^_^,Aij^), {Aij^.Bj^.Ai^)}, A cycle 

(A' j\- j\- \ 
II '2'- ik\^ where the number of 

species nodes defines its order. The set of species nodes in 
a cycle is distinct, but there may be a repetition among the 
reaction nodes. This is because negative paths containing 
the same nodes are considered different depending on the 
starting species node. For example, C = (I3 53) Figure 1 
is a cycle formed by the two negative paths [Ai, ^3, A2] and 

A cycle is positive if it contains an even number of 
negative paths and negative if it contains an odd num- 
ber of negative paths. The sign of a cycle C can also be 
determined by the cycle weight which is a product of all 
corresponding weights of negative and positive paths of C 

Kc= n ^-^jkotji) n ^M'^- (11) 

[Ak,Bj,Ai\eC [Ak,BjAi\eC 

For example, C = (^3'^^) (see Figure 1) is a negative cycle 

of order 2 with weight Kq = -1. The cycle C = i^B^^^^ 
(see Figure 1) is a positive cycle of order 2 with weight 
Kc = 1. 

A subgraph g = {Li, L2, . . . , L^} of G consists of edges 
or cycles L/, / = 1, ... ,5, where each species is the begin- 
ning of only one edge, or one path participating in a cycle. 
In other words, the edges and cycles in a subgraph are 
species mutually disjoint. The number of species nodes in 
a subgraph is defined as its order. The subgraph weight is 
defined using the product of the cycle weights (11) and the 
edges weights (10) of the cycles and edges in^ 

= (-1)' n^^^n^-^^^)' (12) 

Ceg Eeg 



where c is the number of cycles in g. For example, the sub- 
graphs = {[Ai.Bs] , C2 = with weight Kg = -1 is 

shown in Figure 2 (bottom right). 

Since more than one path can exist between species 
nodes via different reaction nodes in a bipartite digraph, 
the number of subgraphs through the same node sets may 
be greater than one. The set of all subgraphs g of order k 
with the same species nodes Vi = {A/^, . . . A/^} and reac- 
tion nodes V2 — {Bj^^ . . . Bjj^} sets is called a fragment 
of order k and is denoted by For a fragment 

^k^h f]) define the number 

as the fragment weight. If /<5^ < 0, then is a critical 
fragment. 

For example, the fragment -Ss (534) is shown in Figure 2 
(top left) together with its three subgraphs ^1 = C3 = 
(lstt)'^2 = {[AMX2 = {ilil)}^ndg, = 
{[Ai,55],[A2,53],[A3,54]}. Each of the first two sub- 
graphs ^1 and g2 contains a positive cycle, and thus 
S3 (534) is a critical fragment since 

Ks, = J2^<s = ^<gi+^<g2+Kg, = -1-1+1 = -1 < 0. 

In [3,22] it is shown that the coefficients of the charac- 
teristic polynomial (8) have the following graph-theoretic 
representation 

E^h • • • l^/V 
Ks,— -> k=h....n, 
. . . Ui, 

(14) 

Note that similar terms in aj^ have been combined using 
summation over the subgraphs of a fragment (13) and (14) 
is in a simplified form. It follows by (14) that the corre- 
spondence between a fragment 5^A:(yJ""y^) ^ non-zero 
term in ak{u,w) is one-to-one. For example, the nega- 
tive coefficient in a^iu, w) given in (9) corresponds to the 
critical fragment ^'3(5 34) shown in Figure 2 (top left). 

Critical fragments corresponding uniquely to negative 
terms in (14) are important for the existence of instabili- 
ties as it is explained next. 

Instability criteria for the Jacobian and the bipartite 
digraph 

Here we summarize classical results from bifurcation 
analysis [1] and more recent results relating graph- 
theoretic methods to instabilities [3,4,22]. 

Multistability often arises from a saddle-node bifurca- 
tion in an ordinary differential equations (ODE) model, 
[1,24]. If a saddle-node bifurcation occurs, then a real 
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Figure 2 Critical fragment and subgraplis of tlie reversible substrate inhibition mechanism. Critical fragment S3 (5 34) and constituent 
subgraplis of tlie reversible substrate inhibition mechanism computed by GraTeLPy. (top left) Critical fragment S3 = (5'34)- (top right) Subgraph 
^3 = {[/\i,55],[/\2,53],[/\3,54]}.(bottom left) Subgraph = C3 = (g^;^^'^^). (bottom right) Subgraph ^2 = {[A],Bs],C2 = (e';^')}- 



eigenvalue k{u,w) of J(u,w) changes sign as the param- 
eters (u, w) change values. Hence, a necessary condition 
for multistability arising from a saddle-node bifurcation is 
ayi(u, w) = det(— /(w, w)) = 0 for some parameter values 
of (u, w) [1]. 

Often ODE models of biochemical mechanisms (4) have 
mass conservation relations reducing the rank r of the sto- 
ichiometric matrix S and the Jacobian J(u, w) to r < n, 
which means that the last non-zero coefficient in (8) is 
ar(u,w). Thus if a saddle-node bifurcation exists, then 
Uriu, w) = 0 for some values of (u, w) [3]. Therefore a crit- 
ical fragment Sr{^j^^""'^0 of order r, corresponding uniquely 
to a negative term in (14) for k = r, is required for a 
saddle-node bifurcation, and thus for multistability [3,22]. 
Thus the potential of a biochemical mechanism (1) for 
multistability depends on the structure of its bipartite 
digraph. 

Oscillations in ODE models of biochemical mechanisms 
(1) often arise from Hopf bifurcation. It is shown in [3], 
that if a coefficient aj^(u,w) > 0, /c g {1, . . . , n — 1} is close 
to zero, then it is possible to choose parameter values for 
(u, w) such that oscillations arising from Hopf bifurcation 
occur. 

The existence of a critical fragment of order 

/c G {1, . . . , w — 1} makes it possible to minimize ai^(u, w) > 
Oj<<n for some parameter values of (u, w) by increasing 



the magnitude of the corresponding negative term in 
UkiUyW), If there are mass conservation relations reduc- 
ing the rank of the Jacobian matrix to r < a critical 
fragment of order /c < r is required to detect 

possible oscillations in an ODE model (4) of a biochem- 
ical mechanism (1). Thus, the existence of oscillations in 
the ODE model of a biochemical mechanism (1) can also 
be determined by the structure of the bipartite digraph. 

Patterns in a corresponding reaction-diffusion model 
to (4) usually arise as a result of Turing instability, Turing 
instability arises when a spatially homogeneous equilib- 
rium is asymptotically stable in the absence of diffusion 
and becomes unstable when diffusion is added to the 
model [25]. For the existence of Turing instability, we 
study the matrix J{Uy w)—iiD, where J(u, w) is the Jacobian 
matrix (6), D is a diagonal matrix with positive diffusion 
coefficients di > 0, i = 1, ... ,n on the diagonal and /x > 0 
is a parameter (/x represents an eigenvalue of the negative 
Laplacian) [25]. Turing instability is associated with a real 
eigenvalue of the matrix J{u, w) — fiD passing through zero 
from left to right as parameter values are varied. In [4,22] 
it is shown that a necessary condition for Turing insta- 
bility is the existence of a critical fragment of 
order k < n. Thus, the potential of a biochemical mecha- 
nism to display Turing instability can be inferred from the 
structure of its bipartite digraph. 
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Implementation 

Recall that the existence of a critical fragment in 
the bipartite digraph of a biochemical mechanism, where 
r is the rank of the stoichiometric matrix Sy can induce 
multistability. Similarly a critical fragment of 
order k < n can induce Turing instability or even oscilla- 
tions. On the other hand if no critical fragments of order 
r are found, then the existence of multistability can be 
excluded for any values of the parameters. Similarly if no 
critical fragments of order k < n are found, then the exis- 
tence of Turing instability can also be excluded for any 
values of the parameters. 

GraTeLPy enumerates all critical fragments of user- 
defined order k for a given biochemical mechanism, thus 
providing the user with information on the potential of a 
biochemical mechanism for multistability, oscillations or 
Turing instability. 

We present in Figure 3 a flowchart that describes 
schematically the algorithm implemented by GraTeLPy. 



First a biochemical mechanism is read from an user- 
provided input text file and its bipartite digraph is gen- 
erated. Then all fragments of an user-defined 
order k are enumerated and placed in a computational 
queue. Each fragment from the queue will be further pro- 
cessed in order to compute its weight (top diamond nodes. 
Figure 3). 

For each fragment in the queue, a linear 

sequence of operations is carried out (central rectangu- 
lar nodes. Figure 3). First all subgraphs ^ of a fragment 
are enumerated, and the weight Kg of each sub- 
graph g is computed. Then the weights of all subgraphs g 
contained in are added to compute the weight 

/<5^ of the given fragment. At this point it is decided based 
on the sign of the weight /<s^, if the fragment is 

critical, i.e., Ks^ < 0 is satisfied. 

Once all of the fragments from the queue have been 
processed, an output based on the potential of the 
biochemical mechanism for some desired instability is 



GraTeLPy Fragment Server 



Read mechanism and generate bipartite digraph. 



Generate all fragments of order k. 



Place all generated fragments in the fragment queue. 



Fragment 1 fetched by a client. Fragment j fetched by a client. Fragment N fetched by a client. 



GraTeLPy Client 



Generate all subgraphs for fragment j . 



Score all subgraphs and produce fragment score. 



Store subgraphs and fragment score received from client. 



Generate output from all results collected. 

GraTeLPy Fragment Server 

Figure 3 Flowchart that summarizes the steps talcen by GraTeLPy to find all critical fragments of a given order. The division of tasl<s 
between a single server and one or more clients is highlighted, (top diamonds) The fragment server reads in the user-specified mechanism file and 
generates the bipartite digraph. The server generates all fragments of an user-defined order k and places them in a queue, (center rectangles) One 
or more client scripts fetch fragments off the queue and process them independently. For each fragment 5^^, a client generates all subgraphs and 
computes the weight of each subgraph. The subgraph weights are then added to compute the weight of the corresponding fragment. The client 
passes the computed data back to the server and fetches another fragment off the queue if the queue is not yet exhausted, (bottom diamonds) 
After preparing the fragment queue, the server waits for the results sent by the clients. Upon receipt of client-computed results for a fragment, the 
server stores these results if the fragment is found to have non-zero weight. Once the queue is exhausted, the server informs the user about the 
number of critical fragments discovered and generates other informative output. 
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created (bottom diamond nodes, Figure 3). The informa- 
tion in the output includes the number of critical frag- 
ments of an user-defined order found by GraTeLPy. Based 
on the number of critical fragments found, GraTeLPy 
states if a biochemical mechanism meets the neces- 
sary condition for multistability or Turing instability, 
and if the mechanism can exhibit oscillations for some 
parameter values. In addition a list of all critical frag- 
ments of a given order detected by GraTeLPy can be 
provided. 

Processing the queue of the enumerated fragments (cen- 
tral rectangular nodes. Figure 3) is inherently parallel as 
each fragment may be handled independently of all other 
fragments. To use this parallelism to our advantage, two 
scripts are created for GraTeLPy that implement a server 
role and a client role, respectively. 

The server script takes care of actions in the top and 
bottom diamond nodes in Figure 3. The server creates 
the bipartite digraph, enumerates all fragments and places 
them in a queue (top diamond nodes). At the end the 
server collects all computed data from the client processes 
before displaying them for the user (bottom diamond 
nodes). 

The client script deals with actions in the central rectan- 
gular nodes in Figure 3. The client fetches a fragment from 
the queue presented by the server, generates all subgraphs 
of the fragment, computes the weights of the subgraphs, 
computes the weight of the fragment, and reports all com- 
puted results back to the server. If the fragment queue 
has not been exhausted, then the client fetches another 
fragment and repeats these steps. 

This server-client architecture allows the user to run 
one or multiple instances of the client script to analyze 
several fragments of a large mechanism in parallel. We 
discuss the technical details of the parallelization in more 
detail in Implementation challenges below. 

In the following subsections we describe in detail the 
implementation of both fragment and subgraph enu- 
meration as these parts presented considerable technical 
challenges during development. 

Fragment enumeration 

It follows by the definition of a fragment given in 
Section The bipartite digraph of a biochemical mecha- 
nism that fragments are identified by the species and reac- 
tion indices of their subgraphs. A fragment ^k{f^'f^""'fj of 
order k contains k unique species indexed by {/i, . . . , 
and k possibly repeated reactions indexed by {/i, ... ,;'/:}. 

Suppose that a given biochemical mechanism has N 
species and R reactions. Using a combinatorial approach, 
we can generate all fragments of order k by pairing the (^) 
unique combinations of species with combinations of 
reaction nodes. This approach generates (^) • R^ possible 



fragments that need to be filtered. This is because many 
of the combinatorially generated fragments do not exist in 
the bipartite digraph of a given biochemical mechanism. 

To save computational time and cost we use a differ- 
ent approach. We note that each fragment ^k{f^'f^""'fj 
contains one subgraph that consists of edges [A/^,^yJ, 
s = 1, . . . , /c. Let us denote by (/ = 1, . . . , m) the 
number of edges that a species A/ in a biochemical mech- 
anism induces, or, is the starting node of. If we assume 
that each species Ai is on average the starting node of 
\E\ = Avg(£/) edges, then this approach generates approx- 
imately (^^) • \E\^ fragments. Empirically, we observe that 
\E\ is usually considerably less than some common values 
for the number of reactions R, Hence this latter approach 
generates fewer fragments than the former combinatorial 
approach. In fact, since fragments correspond uniquely 
to the subgraphs consisting of edges, using this method 
we generate only the fragments that are present in the 
bipartite digraph. 

By using the method of one-to-one correspondence 
between fragments and subgraphs consisting of edges, we 
reduce the number of fragments generated by the com- 
binatorial approach by multiple orders of magnitude. A 
reduction in the number of the generated fragments trans- 
lates directly to a reduction in computational cost. Hence 
the latter approach for fragment generation is an impor- 
tant development in the implementation of GraTeLPy that 
allows for analyzing larger biochemical mechanisms. To 
highlight this reduction in computational cost we plot the 
number of fragments (of varying order) generated with 
both methods for the double-layer mitogen-activated pro- 
tein kinase (MAPK) mechanism in Figure 4. The double- 
layer MAPK mechanism is discussed in more detail in the 
last example in Section Results and discussion. 



E 10^^ 
^ 10^° 










£ 10^ 






1 10' 
^ 10^ 






'3 4 5 6 7 8 9 




order 




Figure 4 Fragment enumeration for double-layer MAPK 




mechanism. Number of fragments of different orders generated for 


the double-layer MAPK network (!) combinatorially (gray) and (ii) 




generated from the unique correspondence between fragments and 


edges-only subgraphs (black). 
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Subgraph enumeration 

Given a fragment ^k(^j\''j^""'j^J we generate all edges 
[Ai^yBj^]y positive paths and negative paths 

[Ai^,Bj^,Ai^], where l,s = 1, . . . , /c, that are induced by the 
species and reactions of the fragment. We will refer col- 
lectively to edges, and positive and negative paths of a 
subgraph as subgraph components. 

The subgraph components of a fragment *^it(yj y'2 y^) 
stored in a lookup table that lists for each species A/^ 
and corresponding reaction Bj^ all subgraph components 
induced by the pair (Ai^.Bj^), The subgraph components 
of a fragment ^k(!^j^f^""^j^^ ^^e generated as follows: 

(i) For a fragment Sk{^-;;;^j^ each edge [A,- ,5yJ 
{s= 1, . . . , A") is identified and stored in the lookup 
table. 

(ii) For each edge [Ai^yBj^] in the lookup table, arcs 
starting at Bj^y such as (Bj^.Ai^) {I = 1, . . . , A") are 
identified. This way all positive paths induced by 
(Ai^.BjJ are generated and added to the lookup table 
as part of the record for species A/^. 

(iii) Similarly to (ii), for each edge [A/^, Bj^] in the lookup 
table, arcs ending at^y^, such as (Ai^^Bp are 
identified. This way all negative paths induced by 
(Ai^,B0 are generated and added to the lookup table 
as part of the record for species Ai^. 

To gain some intuition on how subgraphs can be 
generated, we first describe a simple combinatorial 
approach before we introduce the method implemented 
by GraTeLPy. Suppose that for a fragment ^k{f^'f^""'f^) 
there are {Li^.Li^, . . . subgraph components induced 
by each species {A/^, . . . , A/^}. We can generate combina- 
torially a subgraph^ of ^k{f^'f^""'fj by selecting at random 
one subgraph component per species since each species 
must be the starting node of exactly one component [3]. 
Using this approach we can generate combinatorially all 
possible combinations of subgraph components | • 
• • • \^ik\> ^bat represent all possible subgraph candi- 
dates of a fragment ^k{f^''j^""''j''^)' ^ fragment with \L\ = 
Avg(L/) subgraph components per species on average, this 
method generates \L\^ possible subgraphs. 

Even though this method guarantees that each species 
is the starting node of exactly one subgraph component, 
there may be combinations of paths that do not form 
cycles as defined in Sec. Mathematical background. This is 
because the end species node of a path has to be the start- 
ing species node of another path in a cycle [3] . If we use 
the combinatorial method for generating subgraphs, then 
all candidate subgraphs that do not satisfy the definition of 
a subgraph given in Sec. Mathematical background need 
to be removed which would increase the computational 
cost. 



In the next two subsections we introduce the path graph 
and the cycle graph that will allow us to generate only 
the subgraphs that belong to a given fragment. The imple- 
mentation of the algorithms associated with the path 
graph and the cycle graph by GraTeLPy will allow us to 
further reduce the computational cost. 

Cycle detection: the path graph 

We can avoid generating invalid subgraphs if paths are not 
joined combinatorially, but rather only paths that form 
cycles are joined. Recall that a cycle is a sequence of 
paths where the end species node of each path is the 
starting species node of exactly one other path in the 
sequence. 

Next, we introduce expanded paths, where a negative 
path [Ai, B^,Ay] is converted into two expanded paths 
[A/, B^, Ay] and [Ay, B;^, A/] that are positive. This expan- 
sion is necessary as negative paths can be traversed in both 
directions as explained in Section The bipartite digraph 
of a biochemical mechanism. To enumerate all cycles of 
a given fragment ^k^j^J^""]j^^ we construct the directed 
graph (digraph) O. The nodes of O correspond uniquely 
to the expanded negative paths and the positive paths 
of a fragment Si^f;^';^""'!^). We connect the nodes of the 
digraph O representing paths whose end nodes and start- 
ing nodes are the same. For example, there is a directed 
edge in ^ that starts at a node representing [ A^-^ , By^ , A/2 ] 
and ends at a node representing [Ai2, By2, A/3]. Self-loops 
in ^ from a node back to itself are also permitted and they 
correspond to paths of the form [A/^, By^, A/J. 

To summarize, we generate a digraph O with the follow- 
ing properties: 

• The nodes of O are the expanded negative paths 
and the positive paths of a given fragment 

C (hyi2y-yik\ 

• A directed edge (O/, Oy) exists if and only if the end 
species node of the path corresponding to ^/ is the 
starting species node of the path corresponding to <^j. 
Self-loops (O^-, 00 are permitted and correspond to 
positive paths of the form [Ai^Bj^Ai]. 

We refer to the digraph ^ as the path graph. For a frag- 
ment ^k{f^'^jl""'fj with a total of P paths that include 
all expanded negative paths and all positive paths of 
^^ifih 'fk^' the generation of the path graph O has time 
complexity 0(P(P- 1)). 

To detect the cycles of the path graph O, and ultimately 
the cycles of a given fragment Sk{f^'f^""'f^)y implemen- 
tation of Johnson s algorithm [26] provided by NetworkX 
[27] is used by GraTeLPy. For a fragment ^kif^'jl''''^^) with 
a total number of P expanded negative paths and pos- 
itive paths (corresponding uniquely to the nodes of O), 
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Pe sequential relations between these paths (correspond- 
ing uniquely to the directed edges of O), and Pc cyclic 
sequential relations (corresponding uniquely to the cycles 
of the enumeration of all Pc cycles requires 0((P + 
Pe)(Pc + 1)) units of time [26]. 

Next we illustrate the construction and usage of the path 
graph described above. The path graph <^ for the crit- 
ical fragment S3 (534) (see Figure 2) is shown in Figure 5, 
together with the cycles ci and C2 produced by Johnsons 
algorithm. 

Some of the cycles of O enumerated by NetworkX corre- 
spond to closed paths with revisited nodes in the bipartite 
digraph G, and are therefore not cycles of G. This is the 
case because Johnsons algorithm finds all cycles of all 
lengths of the path graph O. In our current implemen- 
tation, we remove cycles of O that correspond to closed 
paths with revisited nodes of the bipartite digraph G. 
However, further optimization of Johnsons algorithm is 
likely possible, so that only cycles that exist in the bipartite 
digraph G are generated in the first place. 

Cycle combinations: the cycle graph 

Suppose that Pc valid cycles of a given fragment 
^,{[i'''2,--;i/<\ \^2iye been found using the algorithm from the 
previous subsection. Possible candidates for subgraphs of 
*^^0iy'2 a) constructed by creating all possible 

combinations of cycles. In total, there are Xlf^i Ck) P^^' 
sible ways to combine Pc cycles into combinations of k 
cycles with no repeating cycles. Then, edges may have 
to be added to the combinations of cycles in order to 
construct the subgraphs of a fragment ^k{f^'f^""'fj' 

Generally, not all combinations of cycles or edges form 
subgraphs since such combinations may not contain every 
species of a fragment ^k{f^'f^""'fj exactly once. Suppose 
that a given set of cycles has mutually disjoint species 
sets, but the orders of the cycles sum to less than /c. In 




C2 : {[^3,^4, ^2], [^2, 53,^1], [^1,^5, A3]} 

Figure 5 Path graph for the critical fragment of the reversible 
substrate inhibition mechanism. The path graph <I> and the 
detected cycles, that do not have repeated species as starting nodes 
of paths, of the critical fragment S3 = (534) shown in Figure 2. Only 
two cycles ci and cj that are reported previously in [3] are found. 



order to form a subgraph of a fragment ^k{f^'^j^""'fj 
need to amend such a cycle combination with a set of 
edges whose species nodes are in ^k{f^'^jl""'fj> but not in 
any of the cycles. If on average a species A/ is the starting 
node of E edges and if on average we have to add /x edges 
to a cycle combination, then we generate combinatorially 
£^ • Ck) possible subgraphs. 

Many of the combinatorially generated cycle and edges 
combinations will have repeated species nodes, thus ren- 
dering such a combination of edges or cycles invalid as 
a subgraph. Hence we need to verify which of the gen- 
erated combinations of cycles or edges are subgraphs of 
^^^ifi h fk)' validating a subgraph requires 0(l) units 
of time then, on average, validating all possible candidates 
for subgraphs has time complexity 0(£^ • Yl^'Li Ck)) — 
0(Pc-)- reality, validating a combination of cycles or 
edges as a subgraph has greater time complexity than 
0(1). Therefore, the computational cost will be greatly 
reduced if we can generate only subgraphs that require no 
further validation steps. 

We use a similar approach to the one for finding the 
cycles of a given fragment. We will reduce the problem of 
generating cycle or edge combinations forming subgraphs 
to a problem that can be solved with available algorithms 
from the literature. To this end we generate an undirected 
graph r whose nodes correspond uniquely to the cycles of 
a given fragment ^k{f^''jl'''''f^)> that are found using the path 
graph O. Drawing an edge between two nodes (represent- 
ing cycles of ^k{f^'f^""'fj) means that these two cycles do 
not share species nodes and can be combined as a part 
of a subgraph. Next, we formally define the undirected 
graph r 

• A node F/ of F represents a cycle of a given fragment 

o, {h,i2,-,ik\ 
^f<yhh,...,jk)' 

• An edge (F/, Tj) exists if and only if the set of species 
nodes of the cycle represented by F/ and the set of 
species nodes of the cycle represented by Fy are 
disjoint. 

We refer to the undirected graph F defined above as a 
cycle graph. 

If a given set of cycles does not contain a number of 
species equal to the order of the subgraph constructed, 
then species-disjoint edges need to be added. To this end 
the problem of generating a subgraph of a given frag- 
ment ^k{f^'f^""'fj can be reduced to finding all cliques in 
F. Recall that a clique is a set of nodes of an undirected 
graph such that every node is connected to every other 
node from the set [23] . To find all subgraphs of a fragment 
^^(fih fk^' its corresponding undirected graph F should 
be searched for all cliques. This is a standard problem 
in graph theory, known as clique enumeration, that can 
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be solved using existing algorithms from the literature 
[28]. 

As an example, we construct the cycle graph F corre- 
sponding to the fragment ^'3 (534) of the Reversible Inhi- 
bition reaction (2) shown in Figure 2 (top left). We use the 
fact that the cycles ci and C2 of the path graph O, shown in 
Figure 5 have paths that share at least one species. Hence, 
the cycle graph F consists of two nodes corresponding 
to the two valid cycles ci and C2 with no edge connect- 
ing them. Since the cycle graph F constructed from the 
valid cycles in Figure 5 is completely disconnected, we can 
choose one cycle at a time and attempt to construct a valid 
subgraph by adding edges to the cycle. If ci is chosen, 
then the remaining nodes Ai and B5 form the edge [Ai, 
B5] yielding a valid subgraph of order 3, Figure 2 (bottom 
right). If C2 is chosen, then no other nodes remain. Thus 
the cycle C2 forms a valid subgraph of order 3, Figure 2 
(bottom left). 

When generating subgraphs of a fragment combinato- 
rially, the number of subgraph candidates depends on the 
number of subgraph components of a given fragment. 
Using the improved algorithm (based on the path and 
cycle graphs) implemented by GraTeLPy, the number of 
generated subgraphs depends on the number of cycles in 
the path graph. 

To compare the computational cost of the two 
approaches, we count the number of subgraphs generated 
in both cases for 100 randomly selected fragments of vary- 
ing order for the double-layer MAPK network. The results 
are presented in Figure 6, and show that multiple orders of 
magnitude fewer subgraphs are generated by the path and 
cycle graph method in comparison to the combinatorial 
method. 




order 

Figure 6 Subgraph enumeration for double-layer MAPK 
mechanism. Number of subgraph candidates generated for 100 
randomly selected fragments (of indicated order) of the double-layer 
MAPK network. Gray: combinatorial approach, black: using the path 
and cycle graphs. Circles denote averages, squares denote maxima 
(maximal number of subgraph candidates generated for any one of 
1 00 randomly selected fragments). 



Implementation challenges 

As an overarching principle, we have strived to reduce 
code duplication hence we reuse as many components as 
possible from open source libraries. To this end, we have 
used combinatorial standard libraries distributed with the 
programming language Python [29], NetworkX [27] for all 
graph-related operations, and matplotlib [30] for graphi- 
cal output. We note however that matplotlib is an optional 
package and is not required for the core functionality of 
GraTeLPy. 

Over the course of implementation of GraTeLPy we 
encountered combinatorial blowup and memory usage 
as major challenges. Thus we have designed GraTeLPy 
to minimize storage of fragments, subgraphs, and inter- 
mediate structures in memory. To this end we make 
considerable use of the standard Python library itertools 
and the concept of generators that allow us to transfer 
many results from one method to another with minimal 
memory footprint. 

The cycle enumeration method provided by NetworkX 
[27] stores all detected cycles, causing memory shortage 
and overflow due to the large number of generated invalid 
cycles revisiting species. We have amended this library 
method to only store and return valid cycles of the bipar- 
tite digraph, i.e., those cycles that do not revisit species 
nodes. 

Analyzing large networks is computationally unfeasible 
when only a single processor is used. Hence, we have par- 
allelized the code. Pythons global interpreter lock causes 
threaded code to run slowly, so we have used the Mul- 
tiprocessing module [29], which operates by using sub- 
processes rather than threads. We have implemented two 
parallelized versions of the code: 

1. Single multiprocessor machine. Here we use the 
Multiprocessing. Pool system, whereby the 
Multiprocessing module itself launches the requisite 
subprocesses. The code allows the user to specify the 
number of subprocesses that should be run (ideally 
this should match the number of processors available 
on the machine). 

2. Multiple machines client/server. In this case, we 
launch a server process that generates the list of 
fragments to be analyzed. Clients can then be 
launched on any machine with a network connection 
to the server. These clients receive fragments from 
the server and pass back to the server the results of 
the analysis of the fragments. The server collates the 
responses. 

Because the time spent analyzing a fragment is orders 
of magnitude higher than the time required to pass a 
representation of the fragment between a client and the 
server, the parallelization is extremely efficient. We have 
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tested this client/server implementation with over 500 
client processes, and the processing time scales very 
well 

Results and discussion 

GraTeLPy allows the user to enumerate critical fragments 
of an user-defined order. Thus biochemical mechanisms 
can be analyzed for their potential for some instability in 
an efficient way. The existence of multistability requires at 
least one critical fragment of order r, which is the rank of 
the stoichiometric matrix. If a critical fragment of order 
k < n exists, then oscillations may exist for some param- 
eter values. The existence of Turing instability requires at 
least one critical fragment of order k < n, where n is the 
number of species. 

Several examples of biochemical mechanisms of dif- 
ferent sizes are presented in this section. We have 
used GraTeLPy to find the critical fragments of a 
given order in the bipartite digraph of each biochem- 
ical mechanism. The first three examples are smaller 
mechanisms and are used to verify the correctness of 
implementation of GraTeLPy, since their critical frag- 
ments have already been found elsewhere. Furthermore, 
we show that by using GraTeLPy, finding critical frag- 
ments in larger biochemical mechanisms such as the 
MAPK single-layer and MAPK double-layer networks 
becomes feasible. The median running time for finding 
the critical fragments for each biochemical mechanism is 
presented. 

The models and data discussed in this section 
are available at https://github.com/gratelpy/gratelpy- 
supplementary- information. 

Reversible substrate inhibition 

The reversible substrate inhibition model is analyzed for 
multistability in [3] using the graph-theoretic method 
presented here. 

GraTeLPy reads in the biochemical mechanism from a 
plain text file. 

Recall that the bipartite digraph of (15) shown in 
Figure 1. 



The bipartite digraph of (15) contains one critical frag- 
ment 5'3(^'3^) of order 3 found in [3]. GraTeLPy repro- 
duces this fragment and its constituent subgraphs shown 
in Figure 2. 

The median running time with one processor for finding 
the critical fragment ^'3 (534) is 0.05 sec. 

Remark. Note that the critical fragment ^3 = (534) 
corresponds to the negative term in the last non-zero coef- 
ficient (9) of the characteristic polynomial of the Jacobian 
matrix (7). This suggests how we may choose parame- 
ter values for (u, w) so that a saddle-node bifurcation and 
multistability occur. The inequality W4 > wi should be 
satisfied, otherwise a^iu, w) > 0. Also U4 ^ Ui, i = 1, 2, 3 
so that aj,{u, w) is close to zero. In general if is a 

critical fragment, then the species concentrations at equi- 
librium Ui^ > 0, 5 = 1, . . . , /c should be chosen small and 
the rate functions w/^ > 0, 5 = 1, . . . , /c should be chosen 
large in order for a saddle-node bifurcation to occur. 

As a future extension of GraTeLPy, we plan to make 
parameter choices for (u, w) such that some desired insta- 
bility occurs available to the user. 

Glycolysis-Gluconeogenesis switch 

Critical fragments of order smaller than n, the number 
of species in a biochemical mechanism, can induce Tur- 
ing instability in a reaction-diffusion model [4] as well as 
oscillations in an ODE model [3]. 
The biochemical mechanism 
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# Reversible Substrate 

# Inhibition -- Plain 

# text mechanism file 

[Al] -> ; kl 
-> [Al] ; ]<l2 

[Al] + [A2] -> [A3] ; k3 

[A3] -> [A2] ; ]c4 

[Al] + [A3] -> [A4] ; k5 

[A4] -> [Al] + [A3] ; k6 
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represents a glycolysis-gluconeogenesis switch and is 
studied for oscillations in [31]. It has been shown that the 
critical fragments (identified here by GraTeLPy as well) 
are the structural reason for the oscillations [31]. Based 
on the existence of the critical fragments, parameter val- 
ues are chosen such that oscillations occur [31]. Similarly 
the mechanism (15) is studied for the existence of Turing 
instability in [4] . In fact, parameter values are found such 
that Turing instability exists in the reaction-diffusion 
model of (15). 

The bipartite digraph of (15) is shown in Figure 7. 
The stoichiometric matrix associated with the bio- 
chemical mechanism (15) has rank 5. The biochemical 
mechanism meets the necessary criterion for Turing 
instability since critical fragments of order 1 < /c < 5 
exist [4]. In Figure 8 we show the critical fragments 
of order 2 and 3 reported in [4] and identified by 
GraTeLPy. 

The median running time with one processor 
for finding the critical fragments of the bipartite 
digraph of the glycolysis-gluconeogenesis switch is 5.9 
sec. 



Cdc42 network in yeast 

A biochemical mechanism that describes the Cdc42 
dynamics of yeast and cell polarity is studied in [32]. 
The corresponding reaction-diffusion model displays Tur- 
ing instability and patten formation for some parameter 
values [32]. 




Figure 7 Bipartite digraph of the Glycolysis-Gluconeogenesis 
switch mechanism. Bipartite digrapli of tlie glycolysis- 
gluconeogenesis switch. 



The biochemical mechanism of the Cdc42 network is 
given below 



Bi : 
B^: 



E 
E 



RT + En 

1<4 



M, 



M ^ RT + Em. 

B5: En, + RD ^ En,, 

v56 : M + RD ^ M + RT, 



B7: 
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RT % RD, 
Ec + RT ^ M, 



RDL 



k9 



RD + I 



kio 



I + RD, 

RDIm, 



RDIc ^ RDIn,, 

RDIn, ^ RDIn,, 



where RD and RT denote the membrane-bound inactive 
and active form of Cdc42 respectively; I denotes cyto- 
plasmic GDI that forms a membrane-bound complex with 
RD, RDIm, that detaches from the membrane and diffuses 
as RDIc in the cytoplasm. The enzyme E is a complex 
that contains Cdc42-activating Cdc24 and exists in both 
a cytoplasmic and membrane-bound form, Ec and Em, 
respectively. If E is on the membrane, it can form a cat- 
alytic complex M together with RT, that aids activation of 
membrane-bound RD. 

The bipartite digraph of the Cdc42 network is shown 
in Figure 9. The Cdc42 network has a corresponding 
stoichiometric matrix of rank 5. The necessary condi- 
tion for Turing instability requires that a critical fragment 

^^^Ol'yl' a) ^^^^^^^ I < k < 5 exists. GraTeLPy iden- 
tifies 35 critical fragments - among which we find the 
two critical fragments reported in [32] and shown in 
Figure 10. 

The median running time with one processor for finding 
the critical fragments of the bipartite digraph of the Cdc42 
network is 9.7 sec, and with two processors 6.1 sec. 

Single-layer MARK network 

As the size of biochemical mechanisms increases, enu- 
merating the critical fragments of their corresponding 
bipartite digraphs by hand becomes tedious and difficult. 
Using GraTeLPy we can find the critical fragments of a 
given order of larger mechanisms in a short period of time. 
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Figure 8 Critical fragments of the Glycolysis-Gluconeogenesis switch mechanism. Order 2 (bottom) and order 3 (top) critical fragments of tine 
glycolysis-gluconeogenesis switcli, reported in [4] and found by GraTeLPy. 



An example of a larger biochemical mechanism that 
is difficult to analyze by hand is the single-layer MAPK 
network 
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E2, 

Ap + E2, 
ApE2, 
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whose bipartite digraph is shown in Figure 11. The MAPK 
network is a well-known example of a multistable system 
[5,33]. The necessary condition for multstability requires 
the existence of a critical fragment of order equal to the 
rank of the stoichiometric matrix. Since the rank of the 
stoichiometric matrix for the MAPK network equals 6, 
using GraTeLPy, we enumerate all critical fragments of 
order 6. The 9 critical fragments of order 6 of the MAPK 
network are shown in Figure 12. 

The median running time with two processors for find- 
ing the critical fragments of the bipartite digraph of the 
single-layer MAPK network is 10.7 sec. 



Double-layer MAPK network 

For large biochemical mechanisms the number of criti- 
cal fragments of a given order may grow into the dozens 
or hundreds. Thus the task of enumeration by hand of all 
critical fragments of a given order becomes unfeasible, but 
can be accomplished with the help of GraTeLPy. 
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Figure 1 1 Bipartite digraph of the single-layer MARK mechanism. 



An example of a larger biochemical mechanism is pro- 
vided by the double-layer MAPK network which has 12 
species and 18 reactions 
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(16) 



Similarly to the single-layer MAPK network, the double- 
layer MAPK network is known to display multistability. 
The stoichiometric matrix for the double-layer MAPK 
network has rank 9. Therefore the necessary condition 
for multistability requires the existence of at least one 
critical fragment of order 9. GraTeLPy detects 88 crit- 
ical fragments of order 9. The list of all critical frag- 
ments of order 9 of the bipartite digraph of (17) can 
be obtained from https://github.com/gratelpy/gratelpy- 
supplementary- information. 

The median running time of each client for finding the 
critical fragments of the bipartite digraph of the double- 
layer MAPK network is as follows: 141 sec with 100 
clients, 270 sec with 50 clients and roughly 4 hours with a 
single client. 

Conclusions 

We have implemented a graph-theoretic method that 
allows for parameter-free model testing of biochemical 
mechanisms with mass action kinetics for multistabil- 
ity, oscillations and Turing instability. GraTeLPy is open- 
source and is based on a free software. GraTeLPy enables 
users to identify the graph structures referred to as crit- 
ical fragments that can be responsible for the existence 
of some instability in a differential equations model of a 
biochemical mechanism (1). 

At present, GraTeLPy expects that the user converts a 
biochemical mechanism to a text format such as the one 
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Figure 12 Critical fragments of the single-layer MARK mechanism. Critical fragments of tine single-layer MARK network found by GraTeLPy. 
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presented in (15). In a future release we plan to include 
additional functionality so that biochemical mechanisms 
provided in SBML [34] and other formats can be analyzed. 
A list of the critical fragments of a user-defined order is 
provided upon completion. 

We plan a future extension of GraTeLPy where choices 
of parameter values such that some desired instability 
occurs will be offered to the user. This extension will 
be based on the existence of a critical fragment and its 
one-to-one correspondence with a negative term in a 
coefficient of the characteristic polynomial (See Remark 
in the Reversible substrate inhibition Example). 

We also plan to combine GraTeLPy with a new analytic 
method, local perturbation analysis (LPA) [35,36], in order 
to test biochemical mechanisms for pattern formation. 

An extension of GraTeLPy to multigraphs [37] that can 
be used for the analysis of gene regulatory networks [38] 
is also left as a future extension. 

Availability and requirements 

GraTeLPy is available from https://github.com/gratelpy/ 
gratelpy and has the following requirements: Python 2.6 
or 2.7 and NetworkX 1.6 or above. 



Additional file 



Additional file 1 : GraTeLPy Manual: A Practical Software Guide. 
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